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STATEMENT  OF  THE  PROBLEM  STUDIED 

The  research  used  parallelized  computational  fluid  dynamics  methods  to  make  high 
fidelity  simulations  of  thermo-chemical  nonequilibrium  flows  in  rocket  motor  plumes.  This 
work  had  two  main  topics:  recent  thermal  radiative  emission  experimental  data  from  an 
Atlas  launch  were  analyzed  using  a  thermo- chemical  model  with  finite-rate  internal  energy 
relaxation  and  chemical  reactions.  In  addition,  a  detailed  analysis  of  high-temperature  flow 
fields  were  performed  using  a  vibrational  and  electronic  state-specific  excitation  model. 
The  simulation  results  were  then  used  to  compute  the  ultraviolet  and  infrared  radiation 
emitted  from  the  high-temperature  flows  to  compare  with  the  recent  measurements. 


SUMMARY  OF  MOST  IMPORTANT  RESULTS 

This  report  summarizes  the  results  of  this  study,  with  most  focus  on  the  analysis  of 
the  Atlas  II  motor  plumes.  Additional  information  may  be  found  in  the  papers  cited  and 
in  the  recent  Ph.D.  thesis  written  by  Dr.  Krishnendu  Sinha  [1]. 

To  summarize  this  work,  the  most  important  findings  were: 

1.  The  numerical  simulation  of  rocket  motor  plumes  remains  a  very  challenging  task. 
Even  with  advanced  implicit  parallel  methods,  the  combined  effects  of  chemical  reac¬ 
tions,  turbulence,  and  the  vast  range  of  length  scales  makes  the  simulations  difficult. 

2.  The  level  of  turbulent  kinetic  energy  assumed  to  be  in  the  exit  plane  of  the  motor 
exhaust  is  critical  to  predicting  the  level  of  radiative  emission  from  the  plume.  This 
not  so  much  an  effect  of  turbulence  itself,  but  is  due  to  changing  the  total  energy  of 
the  flow  in  the  nozzle  exit  plane. 

3.  For  the  simulations  performed,  other  than  the  effect  just  mentioned,  the  effect  of 
turbulence  on  the  flow  was  fairly  minimal. 

4.  The  numerical  methods  and  gridding  used  for  the  plume  simulations  are  critical.  We 
found  that  third-order  accurate  methods  were  required  to  obtain  sufficiently  accurate 
results,  even  with  multi-million  point  grids.  The  effects  of  the  grid  singularities  on 
the  nozzle  grid  axes  is  a  major  problem  that  has  not  been  overcome  with  existing 
methods.  Solving  this  problem  remains  a  major  challenge  and  should  be  the  subject 
of  a  separate  study. 

5.  The  topology  of  the  rocket  motor  plumes  is  very  dependent  on  the  altitude;  this  effect 
is  reproduced  by  the  simulations,  and  the  change  in  the  nature  of  the  emission  from 
the  flow  fields  is  also  captured. 
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6.  The  low- altitude  (15  km)  simulations  that  we  performed  were  extremely  challenging 
due  to  very  stringent  grid  requirements,  very  thin  shear  layers,  and  rapid  chemical 
kinetics.  The  simulations  resulted  in  plumes  with  a  large-scale  unsteadiness  present. 
This  effect  is  representative  of  the  actual  physical  situation,  though  the  calculations 
were  not  designed  to  capture  unsteady  effects.  See  the  more  detailed  information 
below  for  further  discussion. 

7.  The  presence  of  soot  in  the  motor  exhaust  is  a  critical  element  in  the  generation  of  the 
signature  and  the  change  in  the  signature  from  low  to  high  altitude.  A  more  detailed 
model  of  soot  oxidation  should  be  included  in  the  simulations  to  allow  the  soot  signal 
to  be  reduced  downstream  of  the  nozzle  exit  plane.  This  would  increase  the  accuracy 
of  the  simulations  performed  to  date. 

8.  The  k-e  turbulence  model  equations  are  difficult  to  solve  for  these  high  Mach  number 
flows  with  thin,  high-intensity  shear  layers.  We  modified  the  conventional  linearization 
of  the  source  terms  in  the  k-e  model;  this  resulted  in  at  least  a  factor  of  ten  reduction 
in  the  cost  of  these  calculations. 

9.  The  use  of  the  overlay  method  for  the  simulation  of  complex  chemically  reacting 
flows  is  straight-forward  and  very  powerful.  Essentially  arbitrarily  complex  chemical 
kinetics  models  can  be  solved  to  predict  the  emission  from  trace  species  in  plumes  and 
hypersonic  vehicle  flow  fields. 

The  following  report  discusses  these  findings  in  more  detail.  Additional  information  may 
be  found  in  the  cited  publications. 
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ATLAS  II  ROCKET  MOTOR  PLUME  SIMULATIONS 

Plume  Geometries 

The  plume  geometry  is  the  same  as  that  used  in  previous  papers  [1,2]  The  nozzle  struc¬ 
ture  consists  of  a  linear  arrangement  of  three  nozzles:  a  pair  of  booster  nozzles  surround 
the  low  pressure  sustainer  nozzle.  All  nozzles  have  a  radius  of  0.66  m,  and  are  separated 
by  1.524m  (center-to  center).  Since  all  of  the  nozzles  are  co-linear,  only  one-fourth  of  the 
domain  must  be  modeled  due  to  symmetry.  A  cross-section  of  the  computational  domain 
for  the  Atlas-II  plume  is  shown  in  Fig.  1.  Care  was  taken  to  locate  a  grid  line  at  the  edges 
of  the  booster  and  sustainer  nozzles  in  order  to  reduce  the  numerical  error  at  the  interface 
between  the  nozzle  and  free-stream  flows.  If  the  grid  lines  cross  this  interface,  the  plume 
is  artificially  diffused  and  deformed  due  to  numerical  error. 

The  computational  grid  has  three  grid  blocks:  in  Fig.  1,  the  red  grid  has  54  x  70  x  294 
points,  and  represents  the  sustainer  and  part  of  the  free-stream  flow;  the  green  grid  has 
82  x  22  x  294  points  and  comprises  the  booster  flow;  the  blue  grid  has  102  x  82  x  294  points 
and  represents  the  outer  plume  and  free-stream  flow.  The  latter  grid  is  stretched  away 
from  the  plume  centerline  so  that  the  entire  plume  is  contained  within  the  grid.  The  grid  is 
stretched  in  the  axial  direction  away  from  the  nozzle,  with  some  clustering  in  the  vicinity 
of  the  first  shock  reflection  point  at  the  centerline.  This  4.1  million  point  grid  represents  a 
balance  between  accurately  resolving  the  complex  flow  and  the  computational  cost. 


Symmetry  Plane-2  (m) 

Figure  1 .  Cross-section  of  the  computational  grid  at  the  nozzle  exit  plane. 
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A  single  axisymmetric-equivalent  nozzle  with  the  same  mass  and  energy  flow  is  used 
for  studies  of  grid  resolution,  numerical  methods,  and  turbulence  modeling  parameters. 

Numerical  Method 

The  plume  flow  fields  are  simulated  by  solving  the  three-dimensional  Favre-averaged 
Navier-Stokes  equations  for  reacting  flows  using  the  implicit  full  matrix  Data-Parallel 
Lower-Upper  Relaxation  (DP-LUR)  method  [4],  Finite-rate  chemistry  effects  are  included 
using  a  nine-species  model  (CO2,  H2O,  CO,  N2,  O2,  H2,  OH,  O,  H)  with  ten  reactions  [5]. 

The  simulations  are  for  Atlas-II  flight  conditions  at  two  altitudes:  15  and  40  km.  Free- 
stream,  booster,  and  sustainer  properties  at  the  computed  conditions  are  summarized  in 
Table  1.  The  rotational  modes  are  assumed  to  be  in  equilibrium  with  the  translational 
modes.  Vibrational  nonequilibrium  is  allowed  through  the  solution  of  a  vibrational  energy 
conservation  equation.  Vibrational  relaxation  rates  are  taken  from  Millikan  and  White 
[6]  and  are  used  in  the  standard  Landau-Teller  source  term.  The  Park  two-temperature 
vibration-dissociation  coupling  model  is  used  [7].  Species  viscosities  and  thermal  con¬ 
ductivities  are  calculated  via  least  squares  curve  fits  to  experimental  data  [8,9],  and  the 
mixture  values  are  then  determined  using  Wilke’s  mixing  rule  [10].  Diffusion  coefficients 
are  calculated  using  Fick’s  Law. 


15  km 

40  km 

Single  Nozzle 

Booster 

Sustainer 

Velocity  (m/s) 

535 

1476 

2960 

2930 

3280 

Pressure  (N/m2) 

12100 

278 

68850 

78410 

20200 

Temperature  (K) 

217 

251 

2230 

2380 

1684 

Mass  Fractions: 

C02 

0.000458 

0.000458 

0.302262 

0.324814 

0.337660 

h2o 

0.000000 

0.000000 

0.272443 

0.280699 

0.249950 

CO 

0.000000 

0.000000 

0.412135 

0.381744 

0.396997 

n2 

0.766388 

0.766388 

0.000000 

0.000000 

0.000000 

02 

0.233154 

0.233154 

0.000028 

0.000175 

0.000000 

h2 

0.000000 

0.000000 

0.012194 

0.009945 

0.015381 

OH 

0.000000 

0.000000 

0.000797 

0.002319 

0.000007 

0 

0.000000 

0.000000 

0.000014 

0.000067 

0.000000 

H 

0.000000 

0.000000 

0.000127 

0.000237 

0.000005 

Table  1.  Free-stream  and  nozzle  exit  conditions  at  15  and  40  km. 


The  numerical  fluxes  are  evaluated  using  a  third-order  accurate  MUSCL  [11]  modi¬ 
fied  Steger- Warming  scheme  [12].  Tests  indicate  that  a  second-order  differencing  scheme 
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becomes  unstable  near  the  axis  degeneracies  (located  at  the  nozzle  centerlines  with  the 
present  grid  topology).  This  is  due  to  numerical  error  accumulation  and  amplification 
near  the  axis.  We  have  tested  several  third-order  accurate  differencing  schemes  and  found 
them  to  be  more  suitable  and  stable  near  the  axis.  Hence,  a  third-order  scheme  is  used. 
The  use  of  this  high-order  accurate  differencing  method  significantly  reduces  the  grid  re¬ 
quirements. 

In  all  of  the  cases  presented  in  this  paper  the  numerical  simulation  begins  at  the  exit 
plane  of  the  motor  nozzles,  and  a  simple  supersonic  inflow  boundary  condition  is  assumed. 
A  hyperbolic  tangent  function  is  used  to  smooth  flow  gradients  between  the  nozzle  and 
free-stream.  This  type  of  smoothing  was  previously  shown  to  be  required  to  obtain  a 
grid-independent  solution  [2]. 


Turbulence  Modeling 

The  turbulent  flow  is  simulated  using  the  Favre-averaged  Navier-Stokes  equations 
along  with  the  high  Reynolds  number  version  of  the  k  —  e  turbulence  model  of  Jones 
and  Launder  [13].  The  conservation  equations  for  the  turbulent  kinetic  energy,  k ,  and 
solenoidal  part  of  the  turbulent  dissipation,  es,  are  as  follows. 

dpk  dpujk  d  f  Pt\  dk  __ 

dt  dxj  dxj  \  <Jk )  dxj 

Vk  -  pes  -  ptc  +  p"d" 


dpts 

dt 


dpu3es 

dxj 


d  |Y  pr\  dts_ 

dxj  \  cr€  )  dxj 


cijVk-c2jpe3 

Here,  p  is  the  average  density,  Ui  represents  the  mass-averaged  velocity  and  p,  is  the  molec¬ 
ular  viscosity.  The  turbulent  viscosity,  pr,  is  modeled  as 

pk2 

where  ec  is  the  compressible  dissipation.  Following  the  Sarkar  and  Lakshmanan  model 

[14], 

ec  =  ot\M2e3 


where  M2  =  2k/'yRT  is  the  turbulent  Mach  number.  The  production  of  turbulent  kinetic 
energy,  Vk,  is  given  in  terms  of  the  Reynolds  stress  tensor,  r^, 


Vk 


llj 


dui 

dxj 
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and  the  pressure-dilatation  term,  p"d",  is  modeled  as  [15] 

p"d"  =  -a2M^Vk  +  013  M2pes  . 

In  order  to  improve  predictions  of  axisymmetric  jets,  two  modifications  are  suggested  by 
Dash  [16].  The  first  modifies  the  value  of  the  model  constant  C\,  and  the  second  introduces 
an  extra  term  in  the  e  equation.  Among  the  various  forms  presented  in  literature,  the 
term  proposed  by  Pope  [17]  is  favored  because  of  its  basis  in  vortex  stretching  arguments. 
However,  Lakshmanan  and  Abdol-Hamid  [18]  found  that  these  extra  terms  deteriorate  the 
model  predictions  for  axisymmetric  turbulent  jets.  Therefore,  this  extra  term  is  not  added 
to  the  e-equation  in  the  present  study.  The  model  constants  used  are 

=  0.09,  ci  =  1.60,  C2  =  1.92, 

<7fc  =  1.0,  <r£  =  1.3, 

au  =  1.0,  0:2  =  0.4,  Q!3  =  0.2. 

Axisymmetric  Parametric  Studies 

Turbulent  flow  computations  require  the  free-stream  values  of  k  and  e,  which  in  this 
case  correspond  to  the  exit  plane  of  the  nozzle.  Without  any  experimental  data  for  the 
nozzle  exit  conditions,  we  need  to  choose  reasonable  values  of  the  turbulent  kinetic  energy 
and  dissipation.  In  order  to  assess  the  effect  of  the  inflow  k  and  e  on  the  solution,  we  per¬ 
formed  a  series  of  parametric  studies  using  an  axisymmetric  plume  geometry  at  conditions 
corresponding  the  40  km  altitude. 

First,  we  fixed  the  value  of  k0 0  such  that  the  turbulent  intensity,  T{  —  k/U2,  is  equal 
to  0.01,  and  then  varied  €<*,  over  a  range  of  reasonable  values.  Figure  2  shows  the  contours 
of  the  computed  turbulent  kinetic  energy.  We  see  that  the  flow  with  =  18  x  109  m2/s2 
(which  corresponds  to  pr/ Poo  —  50,  where  poo  is  the  molecular  viscosity  at  the  plume 
exit  plane)  has  very  low  levels  of  k.  However,  the  other  four  turbulent  simulations  have 
significant  levels  of  turbulent  kinetic  energy,  with  most  k  production  occurring  in  the  shear 
layer  and  downstream  of  the  barrel  shock  reflection. 

Comparing  the  five  turbulent  flows,  we  see  that  decreasing  €oo  results  in  increased 
turbulent  kinetic  energy  production,  except  for  the  case  with  the  lowest  level  of  inflow 
turbulent  dissipation.  Thus  there  is  a  threshold  value  of  Coo  required  to  sustain  the  turbu¬ 
lence. 
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0  10 ,  20  30  40 
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Figure  2.  Normalized  turbulent  kinetic  energy  (h/U^)  in  the  axisymmetric  plume  for 
different  values  of  Co©  and  Ti  =  1%. 
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Figure  3.  Temperature  for  the  laminar  flow,  the  laminar  flow  with  modified  plume  inflow 
conditions,  and  a  turbulent  flow  (T,  =  1%,  =  3.0  x  109m2/s2). 

The  temperature  levels  in  the  plume  shear  layer  for  all  the  turbulent  cases  increase 
by  a  few  hundred  degrees  in  comparison  to  the  laminar  plume,  but  there  is  very  little 
variation  between  the  turbulent  cases  themselves.  The  number  density  of  OH  follows  the 
temperature  trend,  thus  we  see  a  higher  OH  concentration  in  the  turbulent  cases,  but 
relatively  little  variation  between  the  turbulent  cases.  From  this  study,  we  conclude  that 
Coo  has  an  effect  on  the  turbulence  level  in  the  plume,  but  the  temperature  and  species 
concentrations  are  not  significantly  affected  by 

In  the  cases  discussed  above,  the  inflow  turbulent  energy  is  taken  to  be  1%  of  the  aver¬ 
age  kinetic  energy  of  the  flow.  Thus,  the  overall  inflow  energy  in  the  turbulent  simulations 
is  higher  than  that  in  the  laminar  computation.  If  we  study  the  turbulent  flows  in  the  re¬ 
gion  close  to  the  inlet,  we  see  that  in  all  cases,  the  turbulent  kinetic  energy  decays  initially 
and  this  extra  energy  is  distributed  to  the  other  energy  modes.  To  assess  the  effect  of  this 
additional  energy,  we  simulated  a  laminar  case  with  the  extra  energy  added  to  the  nozzle 
exit  plane  conditions.  Figure  3  compares  computed  temperature  for  the  regular  laminar 


8 


ARO  Final  Report 


solution,  the  laminar  flow  with  the  additional  inflow  energy,  and  a  turbulent  flow  with 
Coo  =  3.0  x  109m2/s2.  We  see  that  the  increase  in  temperature  in  the  modified  laminar 
case  is  comparable  to  that  in  the  turbulent  case,  showing  that  the  extra  energy  introduced 
in  the  form  of  either  or  modified  nozzle  exit  conditions  is  the  primary  cause  of  the 
increased  temperature  and  the  resulting  higher  product  concentration  in  the  shear  layer. 
Figure  4  plots  the  temperature  and  OH  number  density  across  the  shear  layer  in  these 
three  flows.  Even  a  small  amount  of  extra  energy  over  the  nozzle  exit  conditions  is  enough 
to  cause  a  few  hundred  degrees  temperature  variation  in  the  shear  layer.  As  a  result,  the 
OH  number  density  changes  by  two  orders  of  magnitude.  Therefore,  the  primary  effect  of 
turbulence  in  this  calculation  is  the  added  energy,  rather  than  the  increased  dissipation  in 
the  shear  layer. 

Next,  we  study  the  effect  of  variation  in  k^  on  the  plume  flow  field.  We  find  that 
there  is  a  gradual  increase  in  the  peak  temperature  and  a  corresponding  increase  in  the 
OH  number  density  as  we  increase  koo.  However,  this  increase  is  small  compared  to  that 
shown  in  Fig.  4,  and  thus  increases  in  the  turbulence  intensity  beyond  1%  do  not  have  a 
significant  effect  under  these  conditions.  Finally,  we  study  the  effect  of  grid  spacing  on  the 
k  -  e  turbulence  model.  We  find  that  refining  the  grid  increases  the  level  of  turbulence, 
but  with  a  sufficiently  fine  grid,  the  solution  is  grid  independent. 


Number  density  of  OH  (cm'3) 

10'°  10”  1012  1013  1014  1015 


Figure  4.  Temperature  and  OH  number  density  across  the  shear  layer  in  the  regular 
laminar  flow,  the  laminar  flow  with  modified  inlet  conditions,  and  a  turbulent  turbulent 
flow  (Tt  =  1%,  too  =  3.0  x  109m2/s2). 
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Three-Dimensioned  Atlas-II  Plumes 

In  this  section  we  discuss  the  results  of  numerical  simulations  of  three-dimensional 
Atlas-II  plumes  at  two  flight  conditions:  15  and  40  km  altitude.  The  40  km  results  were 
reported  in  a  previous  paper  [3],  but  it  is  interesting  to  contrast  these  results  with  the 
lower  altitude  results.  Both  flight  conditions  have  been  computed  with  a  laminar  flow 
model  and  the  two-equation  turbulence  model  discussed  above. 

The  results  presented  here  are  primarily  in  the  form  of  temperature  and  OH  number 
density  contours.  This  gives  a  qualitative  description  of  the  flow  fields  and  shows  how 
the  different  flight  conditions  affect  the  overall  flow  topology.  The  results  have  been  used 
to  compute  the  radiative  emission  from  the  Atlas-II  flow  fields  to  compare  with  emission 
measurements  [19]. 

The  turbulent  cases  were  initialized  by  the  turbulence  quantities  that  gave  the  highest 
turbulence  intensities  in  the  shear  layer  for  the  axisymmetric  40  km  plume  (Ti  =  1%  and 
Coo  =  3.0  x  109m2/s2). 

It  is  not  possible  to  perform  a  rigorous  grid  resolution  study  for  this  problem  because 
the  current  4.1  million  point  grid  is  already  expensive  to  use.  This  grid  was  constructed 
based  on  the  results  of  the  grid  convergence  studies  performed  in  the  previous  papers  [2,3] , 
but  with  the  additional  benefit  of  using  a  third-order  accurate  upwind-biased  method, 
rather  than  a  pure  upwind  second-order  method.  Thus,  the  present  results  should  be 
sufficiently  well  resolved. 


15  km  Altitude  Flight  Condition  Results 

Figure  5  plots  contours  of  the  temperature  and  the  logarithm  of  the  OH  number 
density  in  the  two  symmetry  planes;  parts  (a)  and  (b)  show  the  results  in  the  plane  of  the 
nozzle  centerlines;  (c)  and  (d)  in  the  plane  perpendicular  to  the  nozzle  axis  plane. 

Figure  5(a)  shows  that  at  the  15  km  conditions,  the  turbulence  modeling  does  not 
affect  the  spreading  of  the  plume  appreciably.  The  high  temperature  regions  in  the  plume 
are  confined  to  the  shear  layer  between  the  plume  and  the  external  flow,  until  about  70  m 
from  the  nozzle  exit;  at  that  point  the  high  temperature  region  expands  to  almost  the  entire 
plume  core.  This  spreading  of  the  high-temperature  region  appears  to  be  related  to  the 
formation  of  unstable  vortical  structures  downstream  of  the  50  m  point.  It  would  appear 
that  this  instability  should  not  be  permitted  by  the  implicit  time  integration  method. 
However,  we  have  performed  parametric  studies  to  understand  how  the  formation  of  the 
unstable  vortical  structures  is  affected  by  the  time  step  used  in  the  calculation.  We  find 
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that  there  is  no  discernible  difference  between  the  largest  numerically  stable  time  step 
and  significantly  smaller  integration  steps.  This  is  not  surprising  because  the  time  step  is 
governed  by  the  regions  of  the  flow  with  much  finer  mesh  spacing  than  in  the  region  of  the 
instabilities.  Thus,  the  time  integration  is  very  accurate  in  this  region. 

It  is  completely  reasonable  that  the  plume  shear  layers  become  unstable  and  form 
large-scale  structures,  as  seen  in  the  calculations.  However,  these  calculations  were  not 
designed  for  the  simulation  of  shear  layer  instability  growth.  This  would  require  at  a 
minimum  a  full  360°  flow  field  simulation  to  capture  wave  growth,  a  careful  initialization 
of  the  inflow  to  include  fluctuations  at  appropriate  modes,  and  a  low  dissipation  numerical 
method.  Thus  the  present  results  are  of  questionable  validity.  On  the  other  hand,  the 
results  point  to  the  fact  that  the  flow  is  unsteady,  and  needs  to  be  simulated  using  an 
advanced  method  such  as  a  large-eddy  simulation  (LES)  or  a  detached  eddy  simulation 
(DES). 

The  OH  number  density  plot  shows  an  apparently  more  diffuse  plume,  however  this 
is  largely  a  result  of  plotting  the  logarithm  of  noH-  The  substantial  concentrations  of  OH 
are  confined  to  the  plume  core,  and  largely  occur  in  regions  of  elevated  temperature.  The 
presence  of  the  unstable  vortical  structures  is  also  apparent  in  this  plot. 

The  plume  structure  is  substantially  different  in  the  second  symmetry  plane,  as  seen 
in  Figs.  5(a)  and  (b).  The  high-temperature  region  is  more  tightly  confined  to  the  plume 
centerline,  but  the  plume  has  a  larger  spreading  rate  from  the  axis.  This  is  also  reflected 
in  the  OH  concentration,  but  again  plotting  the  logarithm  enhances  the  apparent  plume 
spreading  rate.  The  more  rapid  spreading  is  a  result  of  the  two  high-pressure  booster  motor 
flows  compressing  the  lower  pressure  sustainer  motor  flow;  the  resulting  high  pressure  gas 
expands  off-axis  and  causes  enhanced  plume  growth  in  this  plane. 

The  unstable  nature  of  the  plume  is  again  apparent  in  this  plot,  but  the  shear  layer 
appears  to  be  unstable  closer  to  the  nozzle  exit  plane,  at  about  30  m.  The  large-scale 
high-temperature  off-axis  regions  downstream  of  about  70  m  are  formed  by  the  unstable 
vortical  structures;  these  high- temperature  regions  propagate  downstream  at  about  the 
local  axial  flow  speed. 

These  results  are  reflected  in  the  temperatures  and  OH  number  densities  plotted  in 
cross-sections  of  the  plume,  as  seen  in  Fig.  6.  For  simplicity,  only  the  turbulent  flow  results 
are  plotted  here.  Close  to  the  nozzle  exit  plane  (Fig.  6(a)  at  4.6  m),  the  primary  flow 
structure  are  the  shear  layers  at  the  edge  of  the  booster  motor  flows.  The  sustainer  motor 
flow  is  strongly  compressed  near  the  plume  centerline  at  this  point,  which  is  upstream  of 
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the  primary  shock  reflection  point.  Small-scale  kinks  in  the  shear  layer  can  be  seen  near 
the  interaction  between  the  booster  and  sustainer  motor  flows;  their  cause  is  unknown. 
Further  from  the  nozzle  exit  plane,  the  flow  has  taken  on  a  very  complicated  structure 
with  a  wide  range  of  scales,  as  seen  in  Fig.  6(b)  (located  56  m  from  the  nozzle  axis).  Very 
complicated  shear  layers  have  been  formed  as  the  shocks  reflect  from  the  plume  centerline, 
the  primary  shear  layer,  and  other  flow  structures.  Further  downstream  (Fig.  6(c)),  the 
plume  has  spread  significantly  and  the  range  of  scales  of  the  structures  has  been  reduced 
due  to  mixing  and  diffusion.  As  seen  in  the  streamwise  contour  plots  (Fig.  5),  the  plume 
spreads  more  in  the  transverse  plane,  although  the  higher  temperature  regions  are  more 
aligned  with  the  nozzle  axis  plane. 

40  km  Altitude  Flight  Condition  Results 

Figure  7  plots  the  same  quantities  as  Fig.  5,  but  at  conditions  corresponding  to  40  km 
altitude.  Here  the  external  flow  speed  is  higher,  but  more  importantly,  the  difference 
between  the  nozzle  exit  pressure  and  the  free-stream  is  much  greater.  This  results  in  a 
more  diffuse  shear  layer  and  a  barrel  shock  reflection  point  much  farther  downstream.  The 
plume  is  entirely  steady  for  both  the  laminar  and  turbulent  calculations,  and  the  primary 
difference  between  the  two  calculations  is  probably  due  to  the  additional  inflow  energy 
that  is  represented  by  the  1%  turbulence  intensity.  Again,  we  see  that  the  plume  spreads 
more  rapidly  in  the  off-nozzle-axis  direction;  in  this  case  this  effect  is  amplified. 

Figure  8  plots  the  temperature  and  OH  concentration  in  the  cross-sectional  plane 
located  at  160  m  from  the  nozzle  exit  plane.  We  see  that  the  plume  structure  is  quite 
different  at  the  higher  altitude.  This  figure  clarifies  that  the  plume  spreads  much  more 
extensively  in  the  plane  perpendicular  to  the  nozzle  axis.  There  is  one  large  nearly  circular 
high  temperature  shear  layer  centered  at  the  plume  axis,  and  two  smaller  nearly  circular 
lobes  expanding  in  the  transverse  plane.  Interestingly,  the  peak  temperature  and  OH 
concentration  occur  in  the  secondary  lobes. 
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Figure  5.  15km  altitude  results:  (a)  tei 
axis  plane;  (c)  temperature  and  (d)  OH  co: 
laminar  solution;  bottom:  k  —  e  solution. 
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Figure  6.  Temperature  and  OH  number  density  contours  in  three  cross-sectional  planes 
of  the  Atlas-II  plume  at  15  km:  (a)  4.6  m;  (b)  56  m;  and  (c)  160  m  from  the  nozzle  exit 
plane. 
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Figure  8.  Temperature  and  OH  number  density  contours  in  a  cross-sectional  plane  160  m 
from  the  nozzle  exit  plane  of  the  Atlas-II  plume  at  40  km. 

Summary 

Numerical  simulations  of  Atlas-II  plumes  have  been  performed  at  two  flight  conditions 
using  a  third-order  accurate  upwind-biased  method  with  implicit  time  integration.  A  nine- 
species  finite-rate  chemical  reaction  model  was  used.  These  are  flying  plume  simulations, 
in  that  the  effects  of  the  hardbody  are  not  included.  Care  was  taken  to  align  the  grid  with 
the  nozzle  lip  to  significantly  reduce  the  numerical  spreading  of  the  plume.  This  required 
the  use  of  a  4.1  million  point  grid  composed  of  three  different  grid  blocks.  Even  with  an 
efficiently  parallelized  code,  a  significant  amount  of  computer  time  was  required  to  obtain 
these  solutions.  The  15  km  altitude  case  is  particularly  demanding  because  of  the  thin 
shear  layers,  large  turbulence  model  source  terms,  and  rapid  chemical  reactions  at  these 
conditions. 

Both  laminar  and  two-equation  turbulence  model  simulations  were  performed.  The 
major  difference  between  these  different  transport  models  appears  to  be  the  result  of  adding 
turbulent  kinetic  energy  to  the  plume  inflow  for  the  turbulent  simulations.  The  plumes 
tend  to  spread  more  quickly  in  the  plane  perpendicular  to  the  nozzle  axis  plane;  this  is 
particularly  true  for  the  40  km  case,  where  two  high-temperature  off-axis  lobes  are  formed. 
The  primary  shear  layers  in  the  15  km  altitude  case  become  unstable  in  our  calculations, 
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resulting  in  the  formation  of  large-scale  structures  that  convect  downstream  through  the 
plume  core.  These  structures  are  regions  of  high  temperature  and  OH  concentration, 
which  will  significantly  affect  the  plume  signature.  However,  the  numerical  simulations 
were  designed  to  capture  these  unsteady  effects,  and  an  appropriate  numerical  approach 
should  be  used  for  more  accurate  predictions  of  their  behavior. 

The  results  of  these  simulations  were  used  to  predict  the  emission  from  the  plumes  at 
15  and  40  km  altitudes.  These  results  were  presented  in  [20]  and  a  paper  has  been  prepared 
and  submitted  for  publication  [21]. 
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OVERLAY  APPROACH  FOR  STATE-SPECIFIC  SIMULATIONS 


As  part  of  the  research  project,  an  overlay  method  was  developed  and  tested  for  the 
simulation  of  flows  with  trace  chemical  species.  The  flow  about  a  hypersonic  vehicle  or 
a  rocket  motor  plume  can  be  modeled  in  two  steps.  First,  the  flow  field  is  simulated 
by  solving  the  Navier-Stokes  equations  that  have  been  extended  to  include  the  effects  of 
finite-rate  chemical  reactions  and  internal  energy  relaxation.  This  solution  includes  all  of 
the  chemical  species  that  may  be  present  in  large  amounts,  or  those  that  are  important  in 
the  chemical  kinetics  of  those  major  species.  Then,  the  mass  conservation  equations  for 
the  trace  chemical  species  are  overlaid  on  the  previously-computed  flow  field.  Since  the 
trace  species  are  present  in  very  small  amounts,  they  do  not  affect  the  mass,  momentum, 
or  energy  of  the  bulk  flow. 


The  computational  fluid  dynamics  method  that  was  used  to  simulate  the  flow  of  the 
major  species  in  the  stagnation  region  of  a  re-entry  experiment  and  in  the  plume  of  the 
thrusting  rocket  has  been  discussed  above.  For  example,  consider  a  baseline  model  for  air 
that  is  a  reacting  mixture  of  thermally  perfect  gases  composed  of  the  following  chemical 
species:  N2,  O2,  NO,  N,  and  O.  These  species  are  allowed  to  react  with  one  another  at 
finite  rates.  Here  NO  is  included  as  a  major  species  because  it  may  be  present  in  relatively 
large  amounts  and  because  is  acts  as  a  reaction  intermediary  in  the  Zeldovich  reactions. 
The  internal  energy  is  described  by  translational,  rotational,  and  vibrational  temperatures. 
Finite-rate  relaxation  of  the  rotational  and  vibrational  temperatures  to  the  translational 
temperature  is  also  modeled. 


The  flow  is  assumed  to  be  described  by  the  extended  Navier-Stokes  equations,  which 
for  two-dimensional  flows  may  be  written  in  conservation  law  form  as 


dU  dF  dG 

aF  +  &  +  “  w’ 


(i) 


where  U  is  the  vector  of  conserved  quantities,  F  and  G  are  the  flux  vectors  in  the  x 
and  y  directions,  respectively,  and  W  is  the  source  vector  due  to  chemical  reactions  and 
the  relaxation  of  the  vibrational  and  rotational  temperatures  toward  the  translational 
temperature.  More  details  concerning  the  form  of  the  fluxes  and  the  source  vector  may 
be  found  in  [22].  The  numerical  method  for  the  solution  of  Eq.  (1)  uses  modified  Steger- 
Warming  flux  vector  splitting,  and  the  steady-state  solution  is  obtained  using  implicit 
Gauss-Seidel  line-relaxation  [12].  The  method  has  been  validated  by  comparison  to  a  wide 
range  of  experimental  data.  The  flow  properties  thus  obtained  will  be  hereafter  referred 
to  as  the  “bulk  solution.” 


If  the  argument  is  made  that  the  introduction  of  trace  species  into  the  bulk  flow  does 
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not  affect  the  energy  or  momentum  of  the  bulk  solution,  then  a  simplified  conservation  law 
for  these  “overlay”  species  can  be  written, 


dps 

dt 


d  , 

+  Tx(p’u 


+  PsUs)  +  ^  ( PsV  +  psVs )  = 


Ws 


S  —  1  j  2,  .  .  .  ?  Tig 


(2) 


where  s  is  the  index  of  each  of  the  ns  overlay  species,  ps  is  the  species  mass  density,  ws 
is  the  chemical  source  term,  u  and  v  are  the  mass-averaged  flow  velocities  in  the  x  and 
y  directions,  and  u8  and  vs  are  the  diffusion  velocities  of  species  s.  The  mass-averaged 
velocities  are  obtained  from  the  bulk  solution,  as  are  the  densities  of  the  major  species  and 
the  temperatures  that  appear  in  the  chemical  source  term.  Hence,  once  the  steady-state 
bulk  solution  has  been  obtained,  Eq.  (2)  can  be  solved.  Note  that  the  overlay  approach 
could  also  be  used  to  add  additional  internal  energy  conservation  equations,  provided  that 
the  energy  contained  in  these  modes  is  small  compared  to  the  total  energy. 

To  make  this  decomposition  of  the  flow  field  into  a  bulk  solution  and  an  overlay 
solution,  we  have  assumed  that  the  overlay  species  do  not  alter  the  state  of  the  bulk 
solution  appreciably.  There  are  three  criteria  that  must  be  met  for  this  condition  to  hold. 
First,  the  trace  species  must  be  present  in  small  amounts  so  that  only  a  small  fraction  of 
the  mass  is  missing  from  the  bulk  flow  solution.  Secondly,  the  presence  of  the  trace  species 
cannot  affect  the  diffusion  coefficients  of  the  bulk  species.  A  very  light  overlay  species, 
such  as  H,  could  possibly  alter  the  mass  diffusion  coefficient  of  the  heavy  bulk  species 
if  it  was  present  in  large  mass  fractions.  However,  as  the  overlay  species  are  present  in 
only  trace  amounts,  this  is  not  an  important  effect.  Thirdly,  and  most  importantly,  the 
overlay  species  cannot  affect  the  chemical  kinetics  of  the  bulk  species.  For  example,  NO 
may  be  present  in  only  small  amounts,  but  it  may  be  the  most  efficient  reaction  path  for 
the  dissociation  of  N2.  Therefore,  it  cannot  be  treated  as  an  overlay  species  unless  N2 
dissociation  is  not  important  under  the  conditions  of  the  simulation.  In  practice,  with  the 
appropriate  choice  of  the  overlay  species,  we  can  apply  this  method  when  the  trace  species 
are  present  at  levels  of  less  than  about  10~4  mass  fraction,  which  would  produce  errors  of 
less  than  0.1%  in  the  bulk  solution  properties. 

The  solution  of  Eq.  (2)  is  much  less  costly  than  the  solution  of  Eq.  (1)  would  be  if  both 
the  bulk  and  overlay  species  were  included.  This  is  because  the  cost  of  solving  Eqs.  (1) 
and  (2)  scale  approximately  as  the  square  of  the  number  of  equations  being  solved.  By 
solving  two  smaller  sets  of  equations  the  computational  cost  is  reduced.  Also,  aside  from 
the  source  terms,  Eq.  (2)  is  linear  in  ps  making  it  easier  to  solve.  Also  note  that  the  flux- 
vector  splitting  is  dramatically  simplified  because  the  flux  Jacobian  is  a  diagonal  matrix, 
and  therefore  the  evaluation  of  the  fluxes  is  trivial.  Also,  many  of  the  components  of  the 
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source  terms  can  be  computed  once  and  stored  for  use  in  all  subsequent  time  steps  during 
convergence  to  steady-state.  Therefore,  the  cost  per  iteration  of  the  overlay  equations  is 
considerably  less  than  that  of  the  full  equation  set.  In  practice,  the  overlay  equation  set 
usually  converges  to  a  steady-state  with  about  a  factor  of  ten  fewer  iterations  than  the 
bulk  solution.  Thus,  it  is  possible  to  incorporate  detailed  excited  state  chemistry  models 
in  the  overlay  solution. 

To  illustrate  the  power  of  the  overlay  approach,  the  results  we  summarize  the  results 
of  a  paper  published  under  the  support  of  this  grant  [23].  In  this  work,  the  overlay  method 
was  used  to  model  the  vibrational  state  distributions  of  NO,  CO,  water,  and  CO2,  which  are 
potential  shock-layer  radiators  in  the  midwave  IR.  The  spectral  predictions  show  that  radi¬ 
ation  from  shock  heated  ambient  CO2  will  be  an  important  contribution  at  most  altitudes 
and  speeds  slower  than  3.5km/s.  This  work  included  a  detailed  vibrational  state-specific 
model  of  CO2  excitation.  Comparisons  of  the  spatial  distributions  of  CO2  vibrational 
states  with  a  corresponding  Boltzmann  distribution  at  the  translational  temperature  show 
that  there  are  substantial  differences  in  the  populations.  These  predictions  are  important 
for  the  design  of  an  upcoming  sounding  rocket  experiment. 
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IMPROVEMENTS  TO  TWO-EQUATION  TURBULENCE 
MODELING  FOR  COMPRESSIBLE  FLOWS 

In  order  to  simulate  a  compressible  turbulent  flow  using  the  k  -  e  turbulence  model, 
we  need  to  solve  the  Reynolds-averaged  equations  governing  the  mean  flow  along  with 
the  transport  equations  for  the  turbulent  kinetic  energy,  k,  and  the  dissipation  rate,  e. 
To  account  for  the  turbulence-chemistry  interactions,  we  may  need  to  solve  an  additional 
transport  equation  for  the  energy  variance. 

The  standard  k-e  model,  developed  for  high  Reynolds  number  flows,  was  originally 
used  along  with  wall  functions  to  account  for  the  viscous  effects  near  a  solid  boundary. 
On  the  other  hand,  the  low  Reynolds  number  versions  of  the  model  allow  integration  of 
the  equations  through  the  viscous  near- wall  region,  thus  enabling  computation  of  complex 
flows  in  which  wall  functions  are  not  valid.  However,  a  very  fine  grid  is  needed  near  the  wall 
to  resolve  the  strong  gradients  of  k  and  e  in  this  region,  and  therefore  the  computation 
becomes  expensive.  In  addition,  these  strong  gradients  make  the  non-linear  low  Reynolds 
number  terms  numerically  very  stiff.  As  a  result,  the  time-step  is  limited  to  small  values 
and  the  solution  requires  a  large  number  of  iterations  to  converge.  This  increases  the 
computational  cost  even  further.  As  a  result,  the  simulation  of  turbulent  flows  using  k-e 
models  becomes  very  computationally  intensive. 

There  are  a  wide  variety  of  numerical  methods  used  to  simulate  turbulent  flows  using 
the  k-e  turbulence  model.  The  methods  can  be  broadly  categorized  into  two  groups: 
explicit  methods  and  implicit  methods.  The  explicit  methods  are  simple  but  the  time-step 
is  limited  to  very  small  values.  As  a  result,  the  solution  requires  an  extremely  large  number 
of  iterations  to  reach  the  steady-state.  On  the  other  hand,  the  implicit  methods  circumvent 
the  time-step  restriction  by  linearizing  the  equations.  Larger  time-steps  can,  therefore,  be 
taken  and  the  solution  converges  in  fewer  iterations.  The  convergence  rate  depends  on  an 
accurate  linearization  of  the  equations.  Therefore,  the  convergence  properties  of  a  method 
may  be  degraded  when  the  low  Reynolds  number,  Re ,  terms  are  not  linearized.  Moreover, 
these  terms  are  noted  to  be  numerically  stiff.  Therefore,  not  linearizing  them  is  expected 
to  deteriorate  convergence  significantly. 

Implicit  Formulation 

For  an  explicit  formulation  the  time-step  is  severely  restricted  by  the  CFL  condition 
(At  <  Ax/(\u\  +  a)).  An  implicit  method  circumvents  this  restriction  by  evaluating  the 
fluxes  and  source  terms  at  the  future  time-level  tn+1  instead  of  tn .  The  fluxes  and  the 
source  terms  at  tn+1  are  obtained  by  linearization;  this  process  is  well  documented  in  the 
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literature. 

Here  we  consider  the  source  terms,  which  are  given  by  Pk ,  Ck  and  Ce  .  The  ex¬ 
pressions  for  these  terms  are  quite  extensive  and  linearizing  all  the  terms  is  a  cumbersome 
task.  Moreover,  in  case  of  an  external  flow  with  a  turbulent  boundary  layer,  only  the 
body- normal  derivatives,  d/dr],  are  important.  This  fact  is  used  to  identify  the  dominant 
contributions  to  the  source  terms. 
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Pk 


£‘-2#in5r 


Ce  ca 


2]iptT  ( d2u\ 2 
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We  use  a  Mach  2.244  boundary  layer  as  a  test  case.  Fig.  9  presents  the  convergence 
history  of  the  solution  for  the  four  different  conventional  implicit  methods.  The  conver¬ 
gence  properties  of  Vandromme’s  method  is  very  similar  to  our  method  which  implies  that 
linearizing  ft  does  not  affect  the  convergence  much.  The  variation  of  ]iT  is  important 
for  the  production  of  turbulence,  and  hence  proper  linearization  of  ft  is  expected  to  be 
essential  for  better  convergence.  However,  linearizing  ft  in  Pk  results  in  a  negative  diag¬ 
onal  contribution  to  the  Jacobian  matrix  and  is  destabilizing.  The  reason  for  little  effect 
of  linearizing  ft  on  the  convergence  may  be  because  of  these  two  effects  counteracting 
each  other. 

Comparison  of  the  convergence  histories  of  the  decoupled  method  and  the  current 
method  shows  that  decoupling  the  turbulent  source  terms  from  the  mean  flow  does  not 
affect  the  convergence  rate  significantly.  This  means  that  the  variation  in  the  source  terms 
is  mainly  determined  by  the  development  of  the  turbulence  variables  k  and  e .  However, 
in  case  of  the  method  that  does  not  linearize  Ck  and  Ce ,  convergence  is  much  slower  than 
the  other  three  methods.  This  implies  that  the  low  Re  terms  are  numerically  very  stiff 
and  need  to  be  linearized  in  order  to  achieve  good  convergence  rates. 

In  the  computation  of  the  Mach  2.244  boundary  layer  using  the  conventional  methods 
which  linearize  Ck  and  C€ ,  the  solution  requires  about  4000  iterations  to  converge.  The 
computational  time,  although  not  prohibitive,  is  quite  large  for  a  simple  flow  computation 
with  relatively  small  number  of  points  (100  x  100  grid).  Increasing  the  number  of  points 
in  the  z -direction  to  200  and  400  requires  16,000  and  50,000  iterations,  respectively.  Thus, 
the  computational  time  become  substantial  in  these  cases.  Therefore,  the  convergence 
properties  of  these  conventional  methods  are  far  from  satisfactory.  In  order  to  identify  the 
problem  that  causes  poor  convergence,  we  studied  the  development  of  the  solution  in  terms 
of  the  variation  of  the  residual.  In  Fig.  9,  the  initial  drop  in  the  residual  corresponds  to 
the  development  of  a  laminar  boundary  layer  on  the  plate.  The  velocity  gradients  in  the 
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laminar  boundary  layer  trigger  the  turbulent  production  mechanisms  and  large  quantities 
of  turbulent  kinetic  energy  are  produced.  This  corresponds  to  the  rise  in  the  value  of 
the  residual  around  iteration  600.  Once  sufficient  turbulence  has  developed,  the  residual 
starts  decreasing  again  (around  iteration  1300).  Although  the  decay  of  the  residual  is  fairly 
steady,  the  residual  frequently  increases  in  the  form  of  sharp  peaks.  We  will  refer  to  these 
peaks  as  “spikes”  in  the  residual  history.  The  time-step  is  limited  by  these  spikes  and 
therefore,  the  solution  requires  a  large  number  of  iterations  to  converge.  In  addition,  the 
time-step  needs  to  be  continuously  monitored  and  frequently  reduced  (Table  6.1)  so  that 
the  spikes  do  not  become  large  and  cause  the  solution  to  diverge.  These  spikes  are  present 
in  case  of  all  the  methods  used  and  they  result  in  slower  convergence.  Thus,  the  key  to 
improve  the  convergence  rates  is  to  eliminate  the  spikes. 


Figure  9.  Comparison  of  convergence  characteristics  of  several  conventional  methods  in 
case  of  a  Mach  2.244  turbulent  boundary  layer  computation. 
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Figure  10.  Convergence  history  of  the  Mach  2.244  boundary  layer  solution  computed 
using  the  original  method  (the  box  represents  the  magnified  view  shown  in  Fig.  13). 


Figure  1 1 .  Effect  of  the  modifications  on  the  convergence  characteristics  of  conventional 
methods  in  the  computation  of  the  Mach  2.244  boundary  layer. 
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We  stabilized  the  spikes  in  two  ways:  the  central  difference  in  the  Jacobian  of  Ck  is 
replaced  by  stable  one-sided  differences,  and  a  low  value  of  turbulent  dissipation  is  used  in 
the  free-stream.  The  Mach  2.244  boundary  layer  is  computed  using  this  modified  method 
and  the  convergence  history  is  compared  with  that  of  the  original  method  in  Fig.  10.  We 
can  see  that  the  spikes  in  the  residual  history  are  eliminated  and  the  only  rise  in  the  residual 
corresponds  to  the  initial  buildup  of  turbulence.  The  final  time-step  can  be  increased  by 
two  orders  of  magnitude,  and  the  steady-state  solution  is  reached  in  less  than  400  iterations. 
Thus,  the  computational  time  is  reduced  by  a  factor  of  10.  Applying  the  modifications  to 
Vandromme’s  method  and  the  decoupled  method  shows  similar  improvements  (Fig.  11).  In 
these  cases,  the  number  of  iterations  is  reduced  by  factors  of  5  and  8,  respectively.  However, 
the  method  which  does  not  linearize  Ck  and  Ce  does  not  show  any  improvement.  This 
is  because  of  the  fact  that  the  modifications  in  the  Jacobian  of  Ck  does  not  affect  this 
numerical  method. 

These  modifications,  and  in  particular,  the  change  in  the  differencing  method  for 
Ck  ,  to  the  implicit  formulation  result  in  dramatic  reductions  in  the  computational  time 
required  to  obtain  solutions  of  complex  flows.  Moreover,  the  time  step  does  not  have 
to  be  continually  adjusted  to  control  the  formation  of  the  spikes,  resulting  in  a  much 
less  monitoring  of  the  calculation.  This  too  results  in  a  large  reduction  in  the  difficulty  of 
performing  these  calculations.  This  approach  was  used  in  the  k-e  simulations  of  the  plumes 
discussed  previously  in  this  report.  Without  this  approach,  those  calculations  would  have 
been  dramatically  more  difficult  to  obtain.  Further  information  about  the  modifications 
to  the  implicit  method  can  be  found  in  [1]. 
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